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ABSTRACT 

We show how to estimate the covariance of the power spectrum of a statistically- 
homogeneous and isotropic density field from a single periodic simulation, by applying 
a set of weightings to the density field, and by measuring the scatter in power spectra 
between different weightings. We recommend a specific set of 52 weightings containing 
only combinations of fundamental modes, constructed to yield a minimum variance 
estimate of the covariance of power. Numerical tests reveal that at nonlinear scales 
the variance of power estimated by the weightings method substantially exceeds that 
estimated from a simple ensemble method. We argue that the discrepancy is caused by 
beat-coupling, in which products of closely spaced Fourier modes couple by nonlinear 
gravitational growth to the beat mode between them. Beat-coupling appears whenever 
nonlinear power is measured from Fourier modes with a finite spread of wavevector, 
and is therefore present in the weightings method but not the ensemble method. 
Beat-coupling inevitably affects real galaxy surveys, whose Fourier modes have finite 
width. Surprisingly, the beat-coupling contribution dominates the covariance of power 
at nonlinear scales, so that, counter-intuitively, it is expected that the covariance of 
nonlinear power in galaxy surveys is dominated not by small scale structure, but rather 
by beat-coupling to the largest scales of the survey. 

Key words: large-scale structure of Universe - methods: data analysis 



1 INTRODUCTION 



The last few years have seen the emergence of a Standard 
ACDM Model of cosmology motivated by and consistent 
with a wide range of observations, including the cosmic 
microwave background, distant supernovae, big-bang nu- 
cleosynthesis, large-scale structure, the abundance of rich 
galaxy cluste rs, and local measur ements of the Hubble 
constant (e.g. iTegmark et al.l l2004b). 

The power spectrum of fluctuations (of temperature, 
density, flux, shear, etc.) is the primary statistic used to con- 
strain cosmological paramete rs from observations of the cos- 
mic microwave backgroun d l|Sj)ejgel_e^_alj|2 00o1). of galax- 
ies iTeg mark et all l2004aj: | Cole et alJ 12003 : ISanchez et alJ 
1200a lEisenstein et alJ 120051) of the Lyman alpha fores t 
JSeliak et al.ll2005l : iLidz et alJl2005i IViel fc Haehneldl2005D 
and of w eak g ravitational le nsing jHoekstoa^^Jjj200^ 
iPen et alJkooa iTakada fc JainlEool Isheldon et alJl2004h 

From a cosmological standpoint, the most precious data 
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lie at large, linear scales, where fluctuations preserve the 
imprint of their primordial generation. A generic, albeit 
not universal, prediction of inflation is that primordial 
fluctuations should be Gaussian. At large, linear scales, 
observations are consi stent with fluctuations being Gaussian 
jKomatsu et al.ll2003l) . 

However, much of the observational data, especially 
those involving galaxies, lies in the translinear or nonlinear 
regime. It remains a matter of ongoing research to elucidate 
the extent to which nonlinear data can be used to constrain 
cosmology. 

We recently began jRimes fc Hamilton! l2005ah a pro- 
gram to measure quantitatively, from cosmological simu- 
lations, the Fisher information content of the nonlinear 
matter power spectrum (specifically, in the first instance, 
the information about the initial amplitude of the linear 
power spectrum). For Gaussian fluctuations, the power 
spectrum contains all possible information about cosmo- 
logical parameters. At nonlinear scales, where fluctuations 
are non-Gaussian, it is natural to start by measuring 
information in the power spectrum, although it seems 
likely that additional information resides in the 3-point 
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and higher order correlation fu nctions jTakada fc Jainl2004t 
ISefusatti fc Scoccimarroll2005l) . 

Measuring the Fisher information in the power spec- 
trum involves measuring the covariance matrix of power. For 
Gaussian fluctuations, the expected covariance of estimates 
of power is known analytically, but at nonlinear scales the 
covariance of power must be estimated from simulations. 

A common way to estimate the covariance ma- 
trix of a quantity is to measure its covariance over 
an e n semble of comput e r sim u lations l|Md]csin_J^_WlLitgj 
[ 19991: IScoccimarro et aP Il999l IZhan fe Eisensteinl 12003 
IZhlmetHiuTl27)o"5riRlme^ However, a re- 

liable estimate of covariance can be comp utationally expen- 
sive, r equiring many, perhaps hu ndreds iMeiksin fc White! 
Il999l : iRimes fc Hamilton! 12005a!) of realizations. On the 
other hand it is physically obvious that the fluctuations 
in the values of quantities over the different parts of a 
single simulation must somehow encode the covariance of 
the quantities. If the covariance could be measured from 
single simulations, then it would be possible to measure 
covariance from fewer, and from higher quality, simulations. 
In any case, the ability to measure covariance from a single 
simulation can be useful in identifying simulations whose 
statistical properties are atypical. 

A fundamental difficulty with estimating covariances 
from single simulations in cosmology is that the data are 
cor related o v er all scales, from small to large. As described 
by iKunschl (Il989f) . such correlations invalidate some of 
the "jackknife" and "bootstrap" schemes suggested in the 
literature. In jackknife, variance is inferred from how much 
a quantity varies when some segments of the data are kept, 
and some deleted. Bootstrap is like jackknife, except that 
deleted segments are replaced with other segments. 

As part of the work leading to the present paper, we 
investigated a form of the bootstrap procedure, in which 
we filled each octant of a simulation cube with a block of 
data selected randomly from the cube. Unfortunately, the 
sharp edges of the blocks introduced undesirable small scale 
power, which seemed to compromise the effort to measure 
covariance o f power reliabl y. Such effects can be mitigated 
by tapering jKunschl lT989l. However, it seemed to us that 
bootstrapping, like jackknifing, is a form of re-weighting 
data, and that surely the best way to re-weight data would 
be to apply the most slowly possible varying weightings. For 
a periodic box, such weightings would be comprised of the 
largest scale modes, the fundamentals. 

In the present paper, 321 we consider applying an 
arbitrary weighting to the density of a periodic cosmological 
simulation, and we show how the power spectrum (and 
its covariance, and the covariance of its covariance) of the 
weighted density are related to the true power spectrum 
(and its covariance, and the covariance of its covariance). 
We confirm mathematically the intuitive idea that weighting 
with fundamentals yields the most reliable estimate of 
covariance of power. Multiplying the density in real space 
by some weighting is equivalent to convolving the density in 
Fourier space with the Fourier transform of the weighting. 
This causes the power spectrum (and its covariance, and the 
covariance of its covariance) to be convolved with the Fourier 
transform of the square (and fourth, and eighth powers) of 
the weighting. The convolution does least damage when the 



weighting window is as narrow as possible in Fourier space, 
which means composed of fundamentals. 

In S0we show how to design a best set of weightings, by 
minimizing the expected variance of the resulting estimate 
of covariance of power. These considerations lead us to 
recommend a specific set of 52 weightings, each consisting 
of a combination of fundamental modes. 

This paper should have stopped neatly at this 
point. Unfortunately, numerical simulations, de scribed in 
a companion paper iRimes fc Hamilton! H0Q6), revealed 
an unexpected (one might say insidious), substantial 
discrepancy at nonlinear scales between the variance of 
power estimated by the weightings method and the variance 
of power estimated by the ensemble method. In JI]we argue 
that this discrepancy arises from beat-coupling, a nonlinear 
gravitational coupling to the large-scale beat mode between 
closely spaced nonlinear wavenumbers, when the power 
spectrum is measured from Fourier modes at anything other 
than infinitely sharp sets of wavenumbers. Surprisingly, in 
cosmologically realistic simulations, the covariance of power 
is dominated at nonlinear scales by this beat-coupling to 
large scales. 

We discuss the beat-coupling problem in j^j Beat- 
coupling is relevant to observations because real galaxy 
surveys yield Fourier modes in finite bands of wavenumber 
fc, of width Afc ~ 1/R where R is a chararacteristic linear 
size of the survey. 

Section H3 summarizes the results. 



2 ESTIMATING THE COVARIANCE OF 
POWER FROM AN ENSEMBLE OF 
WEIGHTED DENSITY FIELDS 

The fundamental idea of this paper is to apply an ensemble 
of weightings to a (non-Gaussian, in general) density field, 
and to estimate the covariance of the power spectrum from 
the scatter in power between different weightings. This 
section derives the relation between the power spectrum 
of a weighted density field and the true power spectrum, 
along with its expected covariance, and the covariance of its 
covariance. 

It is shown, equations 1341 . 1401 . and 1441 . that 
the expected ((covariance of) covariance of) shell-averaged 
power of weighted density fields is simply proportional 
to the true ((covariance of) covariance of) shell-averaged 
power, provided that two approximations are made. The 
two approximations are, firstly, that the power spectrum 
and trispectrum are sufficiently slowly varying functions of 
their arguments, equations 1171 and 1371 . and, secondly, 
that power is estimated in sufficiently broad shells in fc- 
space, equation 1391 . The required approximations are most 
accurate if the weightings contain only the largest scale 
Fourier modes, such as the weightings containing only 
fundamental modes proposed in JJ] 

As will be discussed in JIJ the apparently innocent 
assumption, equation 1371 . that the trispectrum is a slowly 
varying function of its arguments, is incorrect, because it 
sets to zero some important beat-coupling contributions. 
However, it is convenient to pretend in this section and 
the next, an d SHU that the assumption 1371 is true, 
and then to consider in Jl] how the results are modified 
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when the beat-coupling contributions to the trispectrum 
are included. Ultimately we find, H4.5I that the weightings 
method remains valid when beat-couplings are included, 
and, H4.6I that the minimum variance weightings derived 
in J|| while no longer exactly minimum variance, should be 
close enough to remain good for practical application. 

This section is necessarily rather technical, because it 
is necessary to distinguish carefully between various flavours 
of power spectrum: estimated versus expected; unweighted 
versus weighted; non- shell- averaged versus shell- averaged. 
Subsections 12.11 to 12.71 present expressions for the various 
power spectra, their covariances, and the covariances of 
their covariances. Subsections 12.81 and 12.91 show how the 
expressions are modified when, as is usually the case, 
deviations in power must be measured relative to an 
estimated rather than an expected value of power. 



2.1 The power spectrum 

Let p(r) denote the density of a statistically homogeneous 
random field at position r in a periodic box. Choose the unit 
of length so that the box has unit side. The density p(r) 
might represent, perhaps, a realization of the nonlinear ly 
evolved distribution of dark matter, or of galaxies. The 
density could be either continuous or discrete (particles). 
Expanded in Fourier modes p(k), the density p(r) is 1 



P(r) = ^2p(k)e 



(1) 



Thanks to periodicity, the sum is over an integral lattice of 
wavenumbers, k — {k x ,k y ,k z } with integer k x , k y , k z . 

The expectation value (p(r)) of the density defines the 
true mean density p, which without loss of generality we 
take to equal unity 

p={p(r)) = l. (2) 
The deviation Ap(r) of the density from the mean is 



Ap(r) = p(r) - p 



(•3) 



The expectation values of the Fourier amplitudes vanish, 
(p(k)) — 0, except for the zero'th mode, whose expectation 
value equals the mean density, (p(0)) = p. The Fourier 
amplitude p(0) of the zero'th mode is the actual density 
of the realization, which could be equal to, or differ slightly 
from, the true mean density p, depending on whether the 
mean density of the realization was constrained to equal the 
true density, or not. 

Because the density field is by assumption statistically 
homogeneous, the expected covariance of Fourier amplitudes 
p(k) is a diagonal matrix 



(Ap(fci) Ap(fc 2 )} = l fel+fe2 P(fci) . 
Here lk denotes the discrete delta-function, 



l fc = 



1 if k = 

otherwise 



(4) 



(5) 



1 The same symbol p is used in both real and Fourier space. 
The justification for this notation is that p is the same vector in 
Hilbert space irrespective of the basis with r espect to which it 
is expanded. See for example Hamilton 1 2005) for a pedagogical 
exposition. 



and P(k) is the power spectrum. Note that there would 
normally be an extra factor of p -2 on the left hand side of 
equation JJJ, but it is fine to omit the factor here because 
the mean density is normalized to unity, equation 
The reason for dropping the factor of ~pT 2 is to maintain 
notational consistency with equation (1151 below for the 
power spectrum of weighted density (where the deviation 
in density is necessarily not divided by the mean). 

The symmetry P(—k) — P(k) in equation @ expresses 
pair exchange symmetry. Below, H2.5I we will assume that 
the density field is statistical isotropic, in which case the 
power is a function P(k) only of the scalar wavenumber k 
k . but for now we stick to the more general case where 
power is a function P(k) of vector wavenumber k. 

2.2 The power spectrum of weighted density 

Let Wi(r) denote the i'th member of a set of real- valued 
weighting functions, and let pi(r) denote the density 
weighted by the i'th weighting 



Pi(r) = Wi(r)p(r) . 



(6) 



The Fourier amplitudes pi{k) of the weighted density are 
convolutions of the Fourier amplitudes of the weighting and 
the density: 



pi{k) = Wi(k')p(k — k') 



h 



Reality of the weighting functions implies 
Wi(-k) = w*(k) . 



(7) 



(8) 



The expected mean p^r) of the weighted density is 
proportional to the weighting, 



PiiT) = (Pi(r)) = Wi(r) 



(9) 



in which a factor of p on the right hand side has been omitted 
because the mean density has been normalized to unity, 
equation The deviation Apj(r) of the weighted density 
from the mean is 



Api(r) ee Pi (r) - Pi (r) 



(10) 



In Fourier space the expected mean p^(fe) of the weighted 
density is 



Pi( k ) = <Pi( fe )) = 



(11) 



and the deviation Ap;(fc) of the weighted density from the 
mean is 



Ap l (fc)EEp l (fc)-p i (fe) 



(12) 



The deviations Api(fc) in the Fourier amplitudes of the 
weighted density are convolutions of the weighting and the 
deviation in the density 



Api(fe) = ^ Wi(fe')Ap(fe - k') 



(13) 



similarly to equation J7J. 

The expected covariance between two weighted densi- 
ties pi(fei) and pj(fc2) at wavenumbers fci and fe 2 is, from 
equations an d 113H . 

(Ap i (fei)Ap J (fc 2 )> =^2w i (k')w j (k 1 +k2-k')P(ki-k') .(14) 
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The weighting breaks statistical homogeneity, so the 
expected covariance matrix of Fourier amplitudes pi(k), 
equation l|14|l . is not diagonal. Nevertheless we define the 
power spectrum Pj(fc) of the i'th weighted density by the 
diagonal elements of the covariance matrix, the variance 



P(fc) = Ap(fc)Ap(-fc) 



(21) 



P(fc) = <A Pl (fc) A Pl (-fc)> 



(15) 



Note that this definition 11511 of the power spectrum Pi(fc) 
differs from the usual definition of power in that the 
deviations Api(k) on the right are Fourier transforms of the 
deviations Api(r) not divided by the mean density Pj(r) = 
Wi(r) (dividing by the mean density would simply unweight 
the weighting, defeating the whole point of the procedure). 
The power spectrum Pi(fc) defined by equation Ijl5^ is 
related to the true power spectrum P(fc) by, equation 114H . 



Pi(fe) = ^2\ Wi (k')\ 2 P(k-k') 



(16) 



Now make the approximation that the power spectrum 
P(fe — fc') at the wavenumber fc — fc' displaced by fc' from fc 
is approximately equal to the power spectrum P(fc) at the 
undisplaced wavenumber fc 

fc') 



P(fc 



P(fc) 



(17) 



This approximation is good provided that the power spec- 
trum P(fc) is slowly varying as a function of wavenumber 
fc, and that the displacement fc' is small compared to fc. 
In g3| we constrain the weightings Wj(fe') to contain only 
fundamental modes, fc' = {k' x , k' y , k' z } with k' x , k' y , k' z — 
0,±1, so that the displacement fc' is as small as it can be 
without being zero, and the approximation 1171 is therefore 
as good as it can be. The approximation I17H becomes exact 
in the case of a constant, or shot noise, power spectrum 
P(fc), except at fc — fc' = 0. 

Under approximation 1171 . the power spectrum of the 
i'th weighted density is 



Pi(fc) wP(fc) V Wfc')| 



(18) 



which is just proportional to the true power spectrum P(fc). 

Without loss of generality, let each weighting Wi(k') 
be normalized so that the factor on the right hand side of 
equation 1181 is unity 



(fc') 



(19) 



Then the power spectrum Pi(fc) of the weighted density is 
approximately equal to the true power spectrum P(fc) 

Pi(fc)«P(fc). (20) 

Thus, in the approximation 1171 and with the normal- 
ization 1191 . measurements of the power spectrum P;(fc) 
of weighted densities provide estimates of the true power 
spectrum P(fc). The plan is to use the scatter in the 
estimates of power over a set of weightings to estimate the 
covariance matrix of power. 

2.3 The covariance of power spectra 

Let P(fc) denote the power spectrum of unweighted density 
at wavevector fc measured from a simulation, the hat 
distinguishing it from the true power spectrum P(fc): 



Below, H2.5I we will invoke statistical isotropy, and we 
will average over a shell in fc-space, but in equation 12H 
there is no averaging because there is just one simulation, 
and just one specific wavenumber fc. Because of statistical 
fluctuations, the estimate P(fc) will in general differ from 
the true power P(fc), but by definition the expectation value 
of the estimate equals the true value, (P(fc)) = P(fc). The 
deviation AP(fc) in the power is the difference between the 
measured and expected value: 



AP(fc) = P(fc) - P(fc) 



(22) 



The expected covariance of power involves the covari- 
ance of the covariance of unweighted densities 

[Ap(fci)Ap(fc 2 ) - l fel+fc2 P(fci)] 

x [Ap(fc 3 )Ap(fc 4 ) P(ks) 

= (Ifc 1 +fe 3 lfe2+fc4 + ^k 1 +k i lk 2 +k 3 ) P{ki)P(k 2 ) 

+ lk 1 +k 2 +k s +k i T(ki,k2, ki, ki) (23) 

which is a sum of a reducible, Gaussian part, the 
terms proportional to P(fci)P(fc2), and an irreducible, 
non-Gaussian part, the term involving the trispectrum 
T(ki, fc2, ka, k4,). Equation 12311 essentially defines what 
is meant by the trispectrum T. Exchange symmetry 
implies that the trispectrum function is invariant under 
permutations of its 4 arguments. The momentum-conserving 
delta-function lfc 1 +fe 2 +fe 3 +fe 4 in front of the trispectrum T 
expresses translation invariance. 

It follows from equation 12311 that the expected 
covariance of estimates of power is 



^AP(fci)AP(fc 2 ) 



(24) 



= (lfc 1+ fe 2 + l fel -fe 2 )P(fci) +T(fc 1 ,-fc 1 ,fc 2 ,-fc 2 ) ■ 

2.4 The covariance of power spectra of weighted 
density 

Similarly to equations 12111 and 12211 . let Pi(fc) denote the 
power spectrum of the i'th weighted density at wavevector 
fc measured from a simulation 



P(fc) = A Pi (k)Api(-k) 



(25) 



and let AP;(fc) denote the deviation between the measured 
and expected value 



APi(fc) = Pi{k)-Pi(k) 



(26) 



The expected covariance between the power spectra of the 
i'th and j'th weighted densities is, from equations 113H and 

lAPiik^APjfa)) = 



w,(fci)wi(fci')wj(fc 2 )wj(fc 2 ) 



fc; +fc' 1 '+fc 2 +fc 2 '=o 

lfci-fc' 1 +fc 2 -fc 2 + lfc 1 -fc' 1 -fc 2 -fc 2 ')P(fci^fci)P( — k\ — k 1 ) 
+ T(fc 1 -fc' 1 ,-fc 1 -fc' 1 ',fc 2 -fc' 2 ,-fc 2 -fc 2 ')l . (27) 



2.5 The covariance of shell-averaged power 
spectra 

Assume now that the unweighted density field ~p(r) i s 
statistically isotropic, so that the true power spectrum P(k) 
is a function only of the absolute value k = |fe| of its 
argument. In estimating the power P(k) from a simulation, 
one would typically average the measured power over a 
spherical shell Vk of wavenumbers in fc-space. Actually the 
arguments below generalize immediately to the case where 
the power is not isotropic, in which case Vk might be chosen 
to be some localized patch in fc-space. However, we shall 
assume isotropy, and refer to Vk as a shell. 

Let p(k) denote the measured power averaged over a 
shell Vk about scalar wavenumber k (the estimated shell- 
averaged power p(k) is written in lower case to distinguish 
it from the estimate P(k) of power at a single specific 
wavevector k): 



(28) 



k£V k 



Here Nk is the number of modes p(k) in the shell Vk- We 
count jo(fe) and its complex conjugate p(—k) as contributing 
two distinct modes, the real and imaginary parts of p(k). 
The expectation value of the estimates p(k) of shell-averaged 
power equals the true shell-averaged power p(k) 



<p(fc)> = P (k) = p ( fe ) 



(29) 



The deviation Ap(k) between the measured and expected 
value of shell-averaged power is 



Ap(k) = p(k) - p(k) = ^~Y1 A ^ 



(30) 



k£V k 



The expected covariance of shell-averaged estimates of 
power is, from equations i-'illt and I2H . 



<Ap(fc 1 )Aj5(fc 2 )) = ^^[2 P ^ 



(31) 



E 



> T(fci,-fci,fe 2 , 

felSV fc , , fc 2 6V t , 



In the usual case, the shells Vk would be taken to be non- 
overlapping, in which case the intersection Vk 1 n Vk 2 in 
equation I31H is equal either to Vu x if V kl and Vk 2 are the 
same shell, or to the empty set if V kl an d Vk 2 are different 
shells. 



2.6 The covariance of shell-averaged power 
spectra of weighted density 

Similarly to equation 1281 . let pi(k) denote the measured 
shell-averaged power spectrum of the i'th weighted density 
at wavenumber k 



Mk) = ^ E Pi W 



(32) 



The expectation value of the estimates Pi(k) is (compare 
eq. ||1|J) 
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In the approximation 11711 of a slowly varying power 
spectrum, and with the normalization 1191 . the expected 
shell-averaged power spectrum Pi(k) of the weighted density 
is approximately equal to the shell-averaged power spectrum 
p(k) of the unweighted density (compare eq. 12011 ) 



Pi(k) fnp(k) 



(34) 



The deviation Api(fc) between the measured and expected 
values is (compare eq. 13011 ) 



Api(k) = pi(k) - pi{k) = — 

JVfc 



feev fc 



AP(k) 



(35) 



The expected covariance of shell-averaged power spectra of 
weighted densities is, from equations I35H and 127L 
i 

(Ap i (k 1 )Apj(k 2 )) 



N kl N k2 
Wi(k'i)wi(k")wj(k' 2 )wj(k2) 



E 



;)P(ki-ki)P(-ki-k'{) 



+ T(fci— ki, —ki—k'l, k2—k' 2 , —k^—k'^) 



(36) 



Assume, analogously to approximation HI 711 for 
the power spectrum, that the trispectrum function 
T(k\ — k'i, — ki — fe", fe 2 — k' 2 , — k-z — k'2) in equation I36H is 
sufficiently slowly varying, and the displacements fej, k", 
k' 2 , k'2 sufficiently small, that 

T(ki — k' 1 , — fei — fe", fc 2 — k' 2 , — fe 2 — k'2) 

w T(fci,-fci,fc 2 ,-fc 2 ) . (37) 

In ||l|we will revisit the approximation 137L and show 
that in fact it is not true, in a way that proves to be 
interesting and observationally relevant. In this section and 
the next, iJ3| however, we will continue to assume that the 
approximation 1371 is valid. 

In the approximations (11711 and 1371 that the power 
spectrum and trispectrum are both approximately constant 
for small displacements of their arguments, the covariance 
of shell- averaged power spectra, equation 136L becomes 
1 

(Api(k 1 )Ap j (k 2 )) 



N kl N k2 



w l (k' 1 )wi(k'i)wj{k'2)Wj(k l 2) 

k' 1 + k '^+ k2 + k ' 2 '=0 

(ifej-fc'j+fca-fe!, + ^k 1 -k' 1 -k 2 -k' 2 ^)P(ki) 

fciev hl , k 2 ev k2 

+ r(fci,-fei,fe 2 ,-fe 2 



(38) 



Consider the Gaussian (P 2 ) part of this expression 138H . 
In the true covariance of shell- averaged power, equation l|31|l . 
the Gaussian part of the covariance is a diagonal matrix, 
with zero covariance between non-overlapping shells. By 
contrast, the Gaussian part of the covariance of power of 
weighted densities, equation l|38fl . is not quite diagonal. In 
effect, the Gaussian variance in each shell is smeared by 
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convolution with the weighting function, causing some of the 
Gaussian variance near the boundaries of adjacent shells to 
leak into covariance between the shells. In we advocate 
restricting the weightings Wi(k) to contain only fundamental 
modes, which keeps smearing to a minimum. Whatever the 
case, if each shell Vk is broad compared the extent of the 
weightings w%(k) in fc-space, then the smearing is relatively 
small, and can be approximated as zero. Mathematically, 
this broad-shell approximation amounts to approximating 

E Ifei-fe'j+fea-fe!, + Ifei-fe'j-feg-*:!,' 

feiev fel , k 2 ev k2 

~ E lfci+fe 2 + lfei-fe 2 = ^ 2. 

fcisv fcl , fc 2 eVfc 2 k 1 ev kl nv k2 

In the broad-shell approximation I39L the expected 
covariance of shell-averaged power spectra of weighted 
densities, equation l|38t . simplifies to 



(Afc(fci)Afe(fe)) « fa (Ap(fci)Aj3(fc 2 )) 
where the factor fa is 



(40) 



fij — 



E 



Wi{k' 1 )wi(k'{)w J (k' 2 )Wj{k2) . (41) 



In real (as opposed to Fourier) space, the factor fa is 
fa = / Wi(r) 2 Wj(rf d 3 r . 



(42) 



Equation l|40|l is the most basic result of the present 
paper. It states that the expected covariance between 
estimates of power from various weightings is proportional to 
the true covariance matrix of power. The nice thing about 
the result (140ft is that the constant of proportionality fa 
depends only on the weightings toj(fe) and Wj(k), and is 
independent both of the power spectrum P(k) and of the 
wavenumbers ki and fe in the covariance (Api(ki)Apj(k2)}- 

2.7 The covariance of the covariance of 

shell-averaged power spectra of weighted 
density 

Equation il I Oil provides the formal mathematical justifica- 
tion for estimating the covariance of power from the scatter 
in estimates of power over an ensemble of weightings of 
density. In J!] we will craft the weightings Wi (k) so as to 
minimize the expected variance of the estimated covariance 
of power. The resulting weightings are "best possible", 
within the framework of the technique. To determine the 
minimum variance estimator, it is necessary to have an 
expression for the (co)variance of the covariance of power, 
which we now derive. 

The expected covariance between estimates 
Api(ki)Api(k 2 ) of covariance of power is a covariance of 
covariance of covariance of densities, an 8-point object. This 
object involves, in addition to the 8-point function, a linear 
combination of products of lower-order functions adding 
to 8 points. The types of terms are (cf. IVerde fc Heavens! 
12001ft 

2 4 , 2-3 2 , 2 2 -4, 2-6, 3-5, 4 2 , 8 (43) 



2 ■ 3 2 signifies a product of a 2-point function with two 
3-point functions, and so on, up to 8, which signifies 
the 8-point function. We do not pause to write out all 
the terms explicitly, because in the same slowly-varying 
and broad-shell approximations that led to equation I40H . 
the covariance of covariance of power spectra of weighted 
densities simplifies to 



[Api(ki)Api(k 2 ) - (A&(&i)Api(fc 2 ))] 

x [Ap J (k 3 )Ap J {k 4 ) - <Ap J (fc 3 )Ap J (fc 4 ))] 
« gy ( [Ap(k 1 )Ap(k 2 ) - (Ap(fci)Ap(fe))] 

x [Ap(k 3 )Ap(k 4 ) - (Ap(k 3 )Ap(k 4 ))] J) (44) 
where gij is, analogously to equation I4H . 



9ij 



E 



(45) 



Wi(k[)wi(k'[)w i (k'2)w i (k'2^Wj(k' s )wj(k'^Wj(k' 4 )wj(k' 4 ) ■ 
In real (as opposed to Fourier) space, the factors gij are 

9a = / Wi(r)V,-(r) 4 d 3 r . (46) 



Equation 14411 states, analogously to equation 14011 . 
that the expected covariance of covariance of power 
spectra of weighted densities is proportional to the true 
covariance of covariance of power. As with the factors 
fa, equation 1411 . the constants of proportionality gij, 
equation 14511 . depend only on the weightings Wi(k) and 
Wj(k), and are independent of the power spectrum P(k) 
or of any of the higher order functions, and are also 
independent of the wavenumbers ki, k 4 in the covariance, 
a gratifyingly simple result. 



2.8 Subtracting the mean power 

The deviation Api(k) of the shell-averaged power spectrum 
of the i'th weighted density was defined above, equa- 
tion 135H . to be the difference between the measured value 
Pi(k) and the expected value Pi(k) of shell-averaged power. 
However, the expected power spectrum Pi(k) (the true 
power spectrum) is probably unknown. Even if the true 
power spectrum is known in the linear regime (because the 
simulation was set up with a known linear power spectrum) , 
the true power spectrum in the non-linear regime is not 
known precisely, but must be estimated from the simulation. 

In practice, therefore, it is necessary to measure the 
deviation in power not from the true value, but rather 
from some estimated mean value. Two strategies naturally 
present themselves. The first strategy is to take the mean 
power spectrum to be the measured power spectrum p(k) 
of the unweighted density of the simulation. In this case 
the deviation Ap[(k) between the measured shell-averaged 
power spectra of the weighted and unweighted densities is 
(the deviation Ap'i(k) is primed to distinguish it from the 
deviation Api(k), eq. 13511 ) 



Ap'i(k) = pi(k) - p(k) 



(47) 



in which 2 signifies a product of four 2-point functions, 



The second strategy is to take the mean power spectrum 
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to be the average over weightings of the measured power 
spectra of weighted densities, iV _1 ^\ pi(fc). In this case 
the deviation Ap i (k) between the measured shell- averaged 
power spectra and their average is (with the same primed 
notation for the deviation ApJ(fc) as in eq. 14711 : it is up to 
the user to decide which strategy to adopt) 



(48) 



The advantage of the first strategy, equation 14711 . is that the 
power spectrum p(k) of the unweighted density is the most 
accurate (by symmetry) estimate of the power spectrum that 
can be measured from a single simulation. Its disadvantage 
is that measurements of power spectra of weighted densities 
yield (slightly) biassed estimates of the power spectrum of 
unweighted density, because the approximation 1171 can lead 
to a slight bias if, as is typical, the power spectrum P(k) is 
not constant. In other words, the approximation Pi(k) ~ 
p(k), equation 1341 . is not an exact equality. Although the 
bias is likely to be small, it contributes systematically to 
estimates of deviations of power, causing the covariance 
of power to be systematically over-estimated. The second 
strategy, equation 1481 . is unaffected by this bias, but 
the statistical uncertainty is slightly larger. Probably the 
sensible thing to do is to apply both strategies, and to check 
that they yield consistent results. 

To allow a concise expression for the covariance of power 
to be written down, it is convenient to introduce Vj(ft), 
defined to be the Fourier transform of the squared real-space 
weighting, Vi(r) = Wi(r) 2 , 



Vi(k) = \_. Wi(k')wi(k" 



(49) 



The normalization condition 1191 on the weightings Wi(k) is 
equivalent to requiring 



Vi(0) = 1 



(50) 



In terms of Vi(k), the factors fij, equation 1411 . relating 
the expected covariance matrix of power spectra of weighted 
densities to the true covariance matrix of power are 



fij =2_^v i (k)v j {-k) 



(51) 



An expression is desired for the covariance of power 
in terms of the deviations Ap'^k), equations 1471 or 1481 . 
instead of Api(k). For this, a modified version of Vi(k) is 
required. For strategy one, equation 1471 . 



(fe = 0) 

Vi(k) (fe £ 0) 



v'i(k) = 

whereas for strategy two, equation 1481 . 

Vi(k) . 



vl(k) = Vi(k) - 



(52) 



(53) 



In either case, the expected covariance ( Ap^(fci) Ap'j (^2)) of 
estimates of shell- averaged power spectra is related to the 
true covariance (Ap(fci) Ap(fc2)} of shell- averaged power by 
(compare eq. 1401 ) 

(A#(fci)A#(fe,)) « f l3 <Ap(fc 1 )Ap(fc 2 )) (54) 
where the factors fL are (compare eq. 1511 1 



fij = y^^(fc)%(~ fe ) 



(55) 



The approximation 1I0 II is valid under the same assumptions 
made in deriving the approximation 1401 . namely the slowly- 
varying approximations 1171 and 1371 . and the broad-shell 
approximation ftfflfl . 



2.9 Subtracting the mean covariance of power 

The expression 1451 for the covariance of covariance of power 
must likewise be modified to allow for the fact that the 
deviations in power must be measured as deviations not from 
the true power spectrum but from either (strategy 1) the 
power spectrum of the unweighted density, or (strategy 2) 
the averaged power spectrum of the weighted densities. 

For this purpose it is convenient to define Ui(k) to be 
the Fourier transform of the fourth power of the real-space 



weighting, m(r) = Vi(r) = u)j(r) 



a(k) 



E 

fe'+fc"=fc 



Vi(k')vi(k") 



(56) 



In terms of Ui(k), the factors gij, equation 1451 . relating 
the expected covariance of covariance of power spectra of 
weighted densities to the true covariance of covariance of 
power are 



9ij 



y\j(fc)Mf(-fc) 



(57) 



To write down an expression for the covariance of the 
covariance of the deviations Ap^(fc) instead of Api(k), define 
a modified version u[(k) of Ui(k) by 



u 'ii k ) = E v 'i( k ') v 'i( k ") 

fc'+fc" = fc 



(58) 



which is the same as equation 1561 but with primed i>;(fc), 
equations 1521 or 1531 . in place of Vi(k). Then the covariance 
of the covariance of the deviations Ap^(fc) is related to the 
true covariance of covariance of shell-averaged power by 
(compare eq. 1441 1 



[A#(fa)A#(*8) - (A$(fci)Ap',(fc 2 )}] 

x [A^(fc 3 )Ap;(fc 4 ) - (A#(ft8)A#(*4)); 



(59) 



(60) 



« g'i^ [Ap(kx)Ap(k 2 ) - (Ap(fc 1 )Ap(fc 2 ))] 
x [Ap(k 3 )Ap{k 4 ) - (Ap(fc 3 )Ap(fc 4 ))] 
where the factors g[j are (compare eq. 1571 1 



fc 



u'i(k)u'j(—k) 



Equation 1591 gives the expected covariance of the dif- 
ference between the estimate of covariance Ap^(fci) Ap^fea) 
and its expectation value (Ap' i (ki)Ap' i (k2)) , but this latter 
expectation value is again an unknown quantity. What 
can actually be measured is the difference between the 
estimate Api(fci)Ap^(fc2) and its average over weightings 
A f ~ 1 y\ Ap' i (ki)Ap' i (k2). To write down an expression for 
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the covariance of the covariance relative to the weightings- 
averaged covariance rather than the expected covariance, 
define a modified version u"(k) of Wj(fe), equation l|58|l . by 

Uj(k) = Uj(k) - u'jjk) . (61) 

i 

Then the covariance of the covariance of the deviations 
Ap'i (k) is related to the true covariance of covariance of shell- 
averaged power by (compare eqs. 14411 and 1591 ') 

( [A#(fci)A#(fc2) - ± Yl ^Pk(ki)Ap' k (k 2 )] 

k 

x [A#(fc 3 )A#(fc 4 ) A PUfc3)Ap;(fc 4 )] ) 

« ( [Ap(fer)Ap(fe) - (Ap(fci)Ap(fe))] 

x [Ap(fc 3 )Ap(fc 4 ) - (Ap(fc 3 )Ap(fe 4 )>] \ (62) 
where the factors g'/j are (compare eqs. (1571 and I6U1 ') 
g" 3 = y^^u"(k)u"(-k) . (63) 

k 

Approximations 1591 and 16211 are valid under the same 
approximations as approximations 1401 and 1441 . namely 
the slowly- varying approximations 1171 and 1371 . and the 
broad-shell approximation 1391 . 



3 MINIMUM VARIANCE WEIGHTINGS 

It was shown in J5| that the expected covariance between 
shell-averaged power spectra of weighted densities is 
proportional to the true covariance of shell-average power, 
equation 1541 . It follows that the scatter in estimates of 
power from different weightings can be used to estimate the 
true covariance of power. In this section we use minimum 
variance arguments to derive a set of 52 weightings, 
equation 1711 . which we recommend, H3.6I for practical 
application. 

In this section as in the previous one, we continue 
to ignore the beat-coupling contributions to the (covariance 
of) covariance of power. These beat-couplings are discussed 
in which in >I4.6I concludes that the minimum variance 
weightings derived in the present section, although no 
longer precisely minimum variance, should be satisfactory 
for practical use. 

3.1 Fundamentals and symmetries 

In the first place, we choose to use weightings lUj(fe) that 
contain only combinations of fundamental modes, that is, 
k = {k x ,k y ,k z } with k x , k y , k z running over 0, ±1. By 
restricting the weightings to fundamental modes only, we en- 
sure that the two approximations required for equation 1541 
to be valid are as good as can be. The first approximation 
was the slowly-varying approximation, that both the power 
spectrum P(k) and the trispectrum T(fei, — fei, fc2, — k 2 ) 
remain approximately constant, equations 1171 and 1371 . 
when their arguments are displaced by the extent of the 
weightings lUj(fe'), that is, by amounts k' for which u>i(k') 
is non-zero. The second approximation was the broad-shell 



approximation, that the shells Vk over which the estimated 
power pi (k) is averaged are broad compared to the extent of 
the weightings Wi(k'), which reduces the relative importance 
of smearing of Gaussian variance from the edges of adjacent 
shells into covariance between the shells. 

In the second place, we choose to use weightings that 
are symmetrically related to each other, which seems a 
natural thing to do given the cubic symmetry of a periodic 
box. Choosing a symmetrically related set of weightings not 
only simplifies practical application of the procedure, but 
also simplifies the mathematics of determining a best set of 
Fourier coefficients Wi(k), as will be seen in H3. 21 below. 

There are 48 rotational and reflectional transformations 
of a cube, corresponding to choosing the x-axis in any 
of 6 directions, then the t/-axis in any of 4 directions 
perpendicular to the :r-axis, and finally the z-axis in either 
of the 2 directions perpendicular to the x- and y-axes. 

To the rotational and reflectional transformations we 
adjoin the possibility of translations by a fraction (half, 
quarter, eighth) of a box along any of the 3 axes, for a 
net total of 48 x 8 3 = 24,576 possible transformations. In 
practice, however, the minimum variance weightings tUj(fe) 
presented in H3.3l prove to possess a high degree of symmetry, 
greatly reducing the number of distinct weightings. 

3.2 How to derive minimum variance weightings 

For brevity, let Xi denote an estimate of the covariance of 
shell-averaged power from the i'th weighted density (the 
arguments fei and k 2 on Xi are suppressed, since they play 
no role in the arguments that follow) 

Xi^A^ffa). (64) 

The quantity /' here is any diagonal element 

/' = fli (65) 

of the matrix of factors /y defined by equation 1551 : the 
diagonal elements fa are identically equal for all i because 
the weightings u>i(k) are by assumption symmetrically 
related. The factor 1/f' in equation 1)640 ensures that Xi 
is, in accordance with equation 1541 , an estimate of the true 
covariance of shell- averaged power, which we abbreviate X, 

{Xi) ^X = (Ap(fci)Ap(fc 2 )) • (66) 

The approximation 1661 is valid under the assumptions 
made in deriving equation 154^ . namely the slowly- 
varying approximations 1171 and 1371 . and the broad-shell 
approximation 1391 . 

Let N denote the number of weightings. Because the 
weightings are by assumption symmetrically related, it 
follows immediately that the best estimate of the true 
covariance of shell-averaged power (Ap(ki) Ap(k2)) will be 
a straight average over the ensemble of weightings 

i 

It remains to determine the best Fourier coefficients 
Wi(k) for a representative weighting i. The best set is that 
which minimizes the expected variance {AX 2 ) = {{X— X)' 2 ) 
of the estimate 1671 . According to equation 1591 . this 
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expected variance (AX 2 ) is approximately proportional to 
a factor that depends on the weightings 

multiplied by another factor that is independent of weight- 
ings, namely the true covariance of covariance of power, the 
expression to the right of the coefficient g'ij in equation 1591 . 
Note that the variance {AX 2 ) is the expected variance 
{{X— X) 2 ) about the true value X, so it is g[j, equation 1601 . 
not g'/j, equation 1631 . that appears in equation 1681 , 

Equation 11681 shows that minimizing the variance 
{AX 2 ) with respect to the coefficients u)j(fe) of the 
weightings is equivalent to minimizing the quantity on 
the right hand side of the proportionality 1681 . From 
equations 1551 . 1581 . and 1601 it follows that this factor can 
be written 

(77^E^ = ^EI u '( fe )l 2 (69) 

ij k 

where u'(fc) denotes the average of u'i(k), equation 1581 . over 
weightings 

= ^E u ^ fe ) • ( 7 °) 

i 

Note that /' = u'(0). Equation l|69|l shows that minimizing 
the variance {AX 2 ) involves computing u'(k), equation 17()|l . 
We evaluate w'(fe) using an algebraic manipulation program 
(Mathematica) as follows. 

A representative weighting Wi(k) contains 27 non-zero 
Fourier coefficients, since by assumption it contains only 
combinations of fundamental modes. The coefficients Wi(k) 
and Wi{—k), which are complex conjugates of each other, 
effectively contribute two coefficients, the real and imaginary 
parts of Wi(k). 

First, evaluate Vi{k), equation 14911 . in terms of the 
coefficients Wi(k) of the representative weighting. The Vi{k) 
are non-zero for 125 values of k, those whose components 
k x , k y , k z run over 0, ±1,±2. Each Vi(k) is a quadratic 
polynomial in the 27 Fourier coefficients. 

Next, modify Vi(k) to get v[(k), equation 1521 . by 
setting the coefficient for k = to zero. Again, each v'^k) 
is a quadratic polynomial in the 27 Fourier coefficients. For 
definiteness, we adopt strategy one, equation 1521 . rather 
than strategy two, equation 1531 . That is, we assume that 
the deviation Ap'^k) in the power spectrum of the i'th 
weighting of density is being measured relative to the power 
spectrum of the unweighted density, rather than relative to 
the average of the power spectra of the weighted densities. 
In the end it turns out, H3.4I that the minimum variance 
solution is the same for both strategies, so there is no loss 
in restricting to strategy one. 

Next, evaluate u'^k), equation 1581 . The u[(k) are non- 
zero for 729 values of k, those whose components k x , k y , k z 
run over 0, ±1, ±4. Each ?4(fc) is a quartic polynomial in 
the 27 Fourier coefficients. 

Next, evaluate u'(k), equation 1701 . the average of ?4(fc) 
over weightings i. Consider first averaging u'i(k) over the 48 
different rotational and reflectional transformations of the 
weighting. The averaged result u'(fc) possesses rotational 
and reflectional symmetry, so that u' (k) is equal to its 



value at k with components permuted and reflected in such 
a way that ^ k x ^ k y ^ k z ^ 4, of which there 
are 35 distinct cases. The rotationally and translationally 
symmetrized function u (k) can be computed by averaging 
the values of u'i{k) at 729 values of k into 35 distinct bins. 

The symmetrized function satisfies u'(fc) = u'(— k), 
so is necessarily real. Thus the absolute value sign around 
it'(fc) 2 in equation 1691 can be omitted. 

Now consider averaging the it; CO over translations by 
half a box in each dimension. There are 2 3 = 8 such 
translations, and each translation is characterized by a triple 
s x , s y , s z giving the number of half boxes translated in 
each dimension, either zero or one for each component. The 
effect of the translation is to multiply each coefficient Wi(k) 
by (-) s * fca:+s !<* ; !' + ^' c *, that is, by ±1 according to whether 
+ s z k z is even or odd. The sign change carries 
through the definitions 1491 of Vi[k) and 1521 of v'i(k) to the 
definition 1581 of u[{k), and thence to the definition 17U1 of 
u (k). That is, the effect of a translation by half a box is to 
multiply «'(fc) by [-y^+s v k y +s z k^ it follows that, after 
averaging over translations, u'(fc) vanishes if any component 
of k is odd, leaving only cases where all components of k 
are even. Consequently, ?4(fe) need be evaluated only at the 
125 wavevectors k all of whose components are even. The 
symmetrized function u'(k) can be computed by averaging 
the values of tt;(fc) at the 125 values of k into the 10 distinct 
bins with even ^ k x ^ k y ^ k z ^ 4. It is amusing 
that increasing the number of weightings (by a factor 8, if 
all translations yield distinct weightings) actually decreases 
the computational work required to find the best Fourier 
coefficients Wi(k). 

Adjoining translations by a quarter of a box simplifies 
the problem of finding the minimum variance solution for 
the coefficients Wi(k) even further. There are 4 3 = 64 such 
translations, and each translation is characterized by a triple 
s x , s y , s z , each component running over to 3, giving 
the number of quarter boxes translated in each dimension. 
The effect of the translation is to multiply each coefficient 
Wi(k) by i s * k x+ a y k y+"z k z . The effect propagates through 
to the symmetrized function u'(k), which is therefore non- 
zero only for the 27 wavevectors k all of whose components 
are multiples of 4. The symmetrized function tt'(fc) can be 
computed by averaging the values of ?4(fc) at the 27 values 
of k into the 4 distinct bins with ^ k x ^ k y ^ k z ^ 4 and 
each component a multiple of 4. 

One more step, adjoining translations by an eighth of a 
box, reduces the problem of finding the minimum variance 
solution to a triviality. After adjoining translations by an 
eighth of a box, the symmetrized function u'(fe) vanishes 
except at k = 0. The function to be minimized, the right 
hand side of equation 1691 . is therefore identically equal to 
1, and any arbitrary weighting therefore yields a minimum 
variance solution. Though amusing, the result is not terribly 
useful, because it involves a vast number, 48 x 8 3 = 24,576, of 
weightings. Physically, if there are enough weightings, then 
together they exhaust the information about the covariance 
of power, however badly crafted the weightings may be. As 
will be seen in H3.3I there are much simpler solutions that 
achieve the absolute minimum possible variance, for which 
the right hand side of equation 1691 equals 1, with far fewer 
weightings. 
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The argument above has shown that the problem of 
finding the minimum variance solution for Wi(k) attains 
its simplest non-trivial form if the weightings are generated 
from a representative weighting by rotations, reflections, and 
translations by quarter of a box, a total of 48 x 4 3 = 3,072 
symmetries. In this case, the weighting-dependent factor in 
the variance of covariance of power, the right hand side of 
equation 16911 . becomes a rational function, a ratio of two 8th 
order polynomials in the 27 Fourier coefficients Wi(k), the 
numerator being a sum ^u'(fc) 2 of squares of 4 quartics, 
and the denominator u (0) 2 the square of a quartic. It is 
this function that we minimize in H3.3l to find a best set of 
weightings. 

The minimum variance solution is independent of the 
overall normalization of the coefficients lOj(fe), since the 
quantity being minimized, the ratio on the right hand side of 
equation 1691 . is independent of the normalization of Wi(k). 
Once the minimum variance solution for the coefficients 
Wi(k) has been found, the coefficients can be renormalized 
to satisfy the normalization condition 1191 that ensures that 
the estimates f>i(k) of the shell-averaged power spectra of 
weighted densities are estimates of the true shell- averaged 
power p(k), equations 1331 and 1341 1. 

3.3 Minimum variance weightings 



The previous subsection, 4'6.2\ described how to obtain the 
coefficients u)j(fe) that minimize the expected variance of the 
estimate of covariance of shell- averaged power that comes 
from averaging over an ensemble of weightings that contain 
only combinations of fundamental modes, and that are 
symmetrically related to each other by rotations, reflections, 
and translations by quarter of a box. 

Numerically, we find not one but three separate sets 
of minimum variance weightings (with hindsight, the sets 
are simple enough that they might perhaps have been 
found without resort to numerics). Each set consists of 
symmetrical transformations of a weighting generated by 
a single mode, namely {1,0,0}, {1,1,0}, and {1,1,1} 
respectively for each of the three sets. Because each 
individual weighting has a rather high degree of symmetry, 
each set has far fewer than the 48 x 4 3 = 3,072 weightings 
expected if all symmetrical transformations yielded distinct 
weightings. Each of the three sets is generated by the 
weighting 



i(fc) = 



±in/a /V2 if k = ±k, 
otherwise 



where ki is one of the three possibilities 



f {1,0,0} 
ki={ {1,1,0} 
I {1,1,1} 



set one: 12 weightings 
set two: 24 weightings 
set three: 16 weightings. 



(71) 



(72) 



In real space, the weighting Wi(r) corresponding to Wi(k) of 
equation I7H is 



Wi(r) = V2c 



2n[ ki.r 



1 

16 



(73) 



The complete set of 12 (24, 16) weightings for each set is 
obtained as follows. In set one (two, three), a factor of 6 
(12, 8) comes from the cubic (dodecahedral, octohedral) 
symmetry of permuting and reflecting the components k x , 




Figure 1. Representative minimum variance weightings Wi(r), 
equation 1731 . for the cases (top) ki = {1,0,0}, and (bottom) 
ki = {1, 1,0}. They are just single Fourier modes, appropriately 
scaled and phased. 



k y , k z of fc, or equivalently the components x, y, z of r. 
A further factor of 2 comes from multiplying Wj(±fc) by ±i, 
equivalent to translating by quarter of a box, or 1/16 — > 5/16 
in equation 1731 . 

The three minimum variance solutions are absolute 
minimum variance, in the sense that each set not only 
minimizes the expression on the right hand side of 
equation 16911 . but it solves u'(k) = for k 7^ 0. This means 
that it is impossible to find better solutions in which all the 
weightings are symmetrically related to each other, which is 
the condition under which equation I69H was derived. 

With the minimum variance solutions in hand, it is 
possible to go back and examine the covariance (Ap'iAp'j), 
equation 154H . between estimates of power from different 
weightings i and j, either within the same set, or across two 
different sets. Estimates of power between two different sets 
are uncorrelated: the covariance (Ap'iAp'j) is zero if i and 
j are drawn from two different sets. If on the other hand 
the weightings i and j are drawn from the same set, then 
it turns out that only half of the weightings, the 6 (12, 8) 
weightings related by the cubic (dodecahedral, octohedral) 
symmetry of permuting and reflecting k x , k y , k z , yield 
distinct estimates of deviation in power. The covariance 
matrix (Ap-Ap^) of estimates of power between the 6 (12, 8) 
cubically (dodecahedrally, octohedrally) related weightings 
is proportional to the unit matrix. However, translating 
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a weighting by quarter of a box, Wj(±k) = ±iWj(±fe), 
yields an estimate of deviation of power that is minus 
that of the original weighting, Ap'j = — Ap[. Actually, this 
is exactly true only if the slowly-varying and thick-shell 
approximations are exactly true (of course, the thick-shell 
approximation is never exactly true). Thus translating a 
weighting by quarter of a box should yield an estimate 
of deviation in power that is highly anti-correlated with 
the original; which should provide a useful check of the 
procedure. 

Translating a weighting by half a box simply changes 
its sign, Wj(±k) = — u>i(±k). This yields an estimate 
of deviation of power that equals exactly (irrespective of 
approximations) that of the original weighting, so yields no 
distinct estimate of deviation in power. These redundant 
translations by half a box have already been omitted from 
the set of 12 (24, 16) weightings. 

The value of /', the factor that converts, equation (1641 . 
estimates Xi of the covariance of power from a weighted 
density field to an estimate of the true covariance of power 



/' = 1/2 



(74) 



the same factor for each of the three sets. 

The expected covariance matrix (AXiAXj) of esti- 
mates of covariance of power equals gL times the true co- 
variance of covariance of power, according to equation 1441 . 
The factors g[j, equation 16011 . are 

C 3/8 if fei = kj 
g^ = l 1/8 iffc i = -fc i (75) 
[ 1/4 otherwise. 

Equation 1751 is valid for weightings i, j both within the 
same set and across different sets. The case ki = kj in 
equation 1751 occurs not only when i = j, but also when 
the weightings i and j are related by translation by quarter 
of a box. The case ki — —kj in equation 1751 occurs not only 
when the weightings i and j are parity conjugates of each 
other, but also when they are parity confugates translated 
by quarter of a box. 

The factors g"j, equation 1631 . which relate the 
covariance {(Xi—X)(Xj—X)) of estimates Xi relative to 
their measured mean X, equation 1671 . as opposed to their 
expected mean X, equation 1661 . are 



1/8 



if k i — k j 



9ij = { -1/8 if k t = -kj 
I otherwise. 



(76) 



An estimate of the uncertainty in the estimate X can 
be deduced by measuring the variance A r_1 '^2 i (Xi— X) 2 in 

the fluctuations about the measured mean X. There is of 
course no point in attempting to estimate the uncertainty 
from N~ 2 Ylij AXiAXj, which is identically zero. The 
true variance {AX 2 } can be estimated from the measured 
variance N' 1 ^(Xi- X) 2 by 



(77) 



N-'E^dij _ 1/4 
N- 1 Ei 9a 1/8 



= 2 



(78) 



in which the factor of 2 comes from (but note the caveat at 
the end of \ 



which corrects for the neglected covariance in the measured 
variance. 



3.4 Minimum variance weightings for strategy two 

The minimum variance weightings derived above assumed, 
for definiteness, strategy one, in which the deviation Ap^(fe) 
in power is taken to be relative to the power spectrum of the 
unweighted density, equation 1471 . An alternative strategy, 
strategy two, is to take the deviation Ap^(fc) in power to be 
relative to the average of the power spectra of the weighted 
densities, equation 1481 . Strategy two yields an estimate of 
covariance of power that has potentially less systematic bias, 
but potentially greater statistical uncertainty. 

As it happens, the minimum variance solution for 
strategy one, H3.3I proves also to solve the minimum 
variance problem for strategy two. Thus the minimum 
variance solution weightings are the same for both strategies. 
Mathematically, expectation values of covariances for the 
two methods differ in that i4(fc) is given for strategy one 
by equation 1521 . and for strategy two by equation 1531 . 
However, for the minimum variance weightings Wi(k) of 
strategy one, equation 1711 and its symmetrical transforma- 
tions, it turns out that N~ £\ v i(k), the term subtracted 
from Vi(k) in strategy two, equation 1531 . is equal to Vj(O) 
if k = 0, and zero otherwise. This is exactly the same as the 
term subtracted from Vi(k) in strategy one, equation 1521 . 
It follows that v'i(k) is the same for the two strategies. 

Although the minimum variance set of weightings is the 
same for both strategies, the two strategies will in general 
yield different estimates of the covariance of power. 



3.5 More minimum variance weightings 

The three minimum variance sets of weightings found 
(numerically) in H3.3l all take the same form, equation 1711 . 
differing only in that they are generated by a different single 
mode, with wavevectors {1,0,0}, {1,1,0}, and {1,1,1} 
respectively. One can check that the result generalizes to 
higher order weightings, in which the wavevector ki in 
equation 1711 is any wavevector with integral components 
(such as {2,0,0}, {2,1,0}, and so on). That is, for any 
wavevector ki with integral components, the weightings 
generated from the weighting of equation 1711 by rotations, 
relections, and translations by quarter of a box, form 
a minimum variance set. All the results of H3.3I (and 
>I3.41 carry through essentially unchanged. In particular, all 
equations 1731 - 1781 remain the same. 

The disadvantage of including higher order weightings 
is that the estimates Xi of the covariance of power become 
increasingly inaccurate as the wavenumber \ki\ of the 
weighting increases, because the slowly- varying approxima- 
tions 1171 and 1371 . and the broad-shell approximation 1391 . 
become increasingly poor as |fej| increases. 

The advantage of including higher order weightings 
is that the more weightings, the better the statistical 
estimate, at least in principle. However, the gain from more 
weightings is no t as great as one m ight hope. Th e Cramer- 
Rao inequality llKendall fc Stuar3ll967t see e.g. iHamilto'nl 
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2005 for a pedagogical derivation) states that the inverse 
variance of the best possible unbiassed estimate X of the 
parameter X must be less than or eq ual to the Fisher 
information F fsee lTeemark et al.lll997l ) in the parameter 
X 



(AX 2 )' 1 sc F 



d 2 ln£ 



(79) 



dX 2 

where C is the likelihood function. To the extent that 
the estimates Xi are Gaussianly distributed (that is, the 
likelihood function is a Gaussian in the estimates Xi 
1 
2 



C 



exp 



-~ ^(AX.AX,)- 1 ^ - X){X } - X) 



(80) 



with covariance (AXiAXj) independent of X), which could 
be a rather poor approximation, the Fisher information F 
in the parameter X approximates the sum of the elements 
of the inverse covariance matrix, 



F « ^ (A^A^) " 



(81) 



In the present case, the covariance matrix (AXiAXj) 
is proportional to <?y, so in approximation that Xi 
are Gaussianly distributed, the Fisher information F is 
proportional to 



p (X 



(82) 



With the coefficients g'ij given by equation 1751 , the quantity 
on the right hand side of equation 1821 proves to be a 
constant, independent of the number of estimates Xi 



(83) 



This constancy of the Fisher information F with respect 
to the number of estimates suggests that there is no gain 
at all in adjoining more and more estimates. However, this 
conclusion is true only to the extent, firstly, that the slowly- 
varying and broad-shell approximations are good, and, 
secondly, that the estimates Xi are Gaussianly distributed, 
neither of which assumptions necessarily holds. All one can 
really conclude is that the gain in statistical accuracy from 
including more estimates is likely to be limited. 

There is however another important consideration 
besides the accuracy of the estimate of the covariance matrix 
of power: it is desirable that the estimated covariance matrix 
be, like the true covariance matrix, strictly positive definite, 
that is, i t should have no zero ( or negative) eigenvalues. As 
noted bv lPan fc Szapudil i2005l) . if a matrix is estimated as 
an average over TV estimates, then its rank can be no greater 
than N. Thus, to obtain a positive definite covariance matrix 
of power for N shells of wavevector, at least N distinct 
estimates Xi are required. 

In j3.6l below we recommend estimating the covariance 
of power from an ensemble of 12 + 24 + 16 = 52 weightings. 
This will yield a positive definite covariance matrix only if 
the covariance of power is estimated over no more than 52 
shells of wavenumber. Since, as noted in H3.3I weightings 
related by translation by quarter of a box yield highly 
anti-correlated estimates of power, hence highly correlated 
estimates of covariance of power, a more conservative 



approach would be to consider that the 52 weightings yield 
only 26 effectively distinct estimates of covariance of power, 
so that the covariance of power can be estimated over no 
more than 26 shells of wavenumber. If (strategy two) the 
deviation of power is measured relative to the measured 
mean over symmetrically related weightings, a (slightly) 
different mean for each of the 3 sets of weightings, then 
3 degrees of freedom are lost, and the covariance of power 
can be estimated over no more than 52 — 3 = 49 shells 
of wavenumber, or more conservatively over no more than 
26 — 3 = 23 shells of power. 

3.6 Recommended strategy 

Here is a step-by-step recipe for applying the weightings 
method to estimate the covariance of power from a periodic 
simulation. 

(i) Select the weightings uti. We recommend the minimum 
variance sets of weightings given by equation 1711 and its 
symmetrical transformations. If the weightings are restricted 
to contain only combinations of fundamental modes, then 
there are three such sets of weightings, equation 172H . and 
the three sets together provide N = 12 + 24 + 16 = 52 
distinct weightings. 

(ii) For each weighting, measure the shell-averaged power 
spectrum Pi(k) of the weighted density field, equations 1321 
and f23t . 

(iii) For each weighting, evaluate the deviation Api(k) 
in the shell-averaged power as the difference between pi(k) 
and, either (strategy one) the shell- averaged power p(k) 
of the unweighted density, or (strategy two) the mean 
N~ 1 '^ / .pi(k) over symmetrically related weightings. The 
advantage of strategy one is that the statistical error is 
potentially smaller, whereas the advantage of strategy two 
is that the systematic bias is potentially smaller. In strategy 
two, it makes sense to subtract the mean separately for 
each symmetrically related set of weightings, because the 
systematic bias is (slightly) different for each set. We 
recommend trying both strategies one and two, and checking 
that they yield consistent results. 

(iv) Estimate the covariance matrix of shell-averaged 
power from the average over all N (52) weightings 



(Ap(ki)Ap{k 2 )) 



N 



^Api(fci)Api(fc2 



(84) 



The factor of 2 in equation 1841 1 is 1//' = 2, equation I74H . 
necessary to convert the average over weightings to an 
estimate of the true covariance of power, equation l|64^ . 



4 BEAT-COUPLING 

This paper should have ended at this point. Unfortunately, 
numerical tests, described in detail in the companion paper 
iRimes fc Hamiltonl|20Q6 }) revealed a serious problem. 

Figure |2] shows the problem. It shows the median 
and quartiles of variance of power measured by the 
weightings method in each of 25 ART ACDM simulations 
of 128/i _1 Mpc box size, compared to the variance of power 
measured over the ensemble of the same 25 simulations. 
Although the two methods agree at linear scales, the 
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Figure 2. Comparison between the normalized variance of 
power measured from 25 ART ACDM simulations by (symbols 
with error bars, indicating median and quartiles) the weightings 
method, and (plain symbols) the ensemble method. The two 
methods disagree substantially at nonlinear scales. Lines show 
the normalized variance predicted by perturbation theory both 
with (solid line), and without (dashed line) the large-scale 
beat-coupling contribution. The dotted line shows the expected 
Gaussian contributio n to the variance. This figu re is a condensed 
version of Figure 5 of lRimes fc Hamilton! <200d) . 

weightings method gives a systematically larger variance at 
nonlinear scales. The discrepancy reaches almost an order of 
magnitude at the smallest s cales measured, fc ~ 5 h~ M pc. 
The reader is referred to iRimes fc Hamilton! (l2006t) for 
details of the simulations and their results. 

This section diagnoses and addresses the problem. The 
next section, discusses the problem and its relevance to 
observations. 



4.1 The cause of the problem: beat-coupling 

The physical cause of the problem illustrated in Figure [5] 
traces to a nonlinear coupling of products of Fourier modes 
closely spaced in wavenumber to the large-scale beat mode 
between them. This beat-coupling, as we refer to it, occurs 
only when power is measured from Fourier modes with a 
finite spread in wavevector, and therefore appears in the 
weightings method (and in observations - see H5.ll below) 
but not in the ensemble method. The beat-coupling is 
surprisingly large, to the point that, as seen in Figure it 
actually dominates the variance of power at nonlinear scales. 

More specifically, in the ensemble method, the power 
spectrum of a periodic simulation is measured from the 
variance Ap(k)Ap(~ fc) of Fourier modes. In the weightings 
method on the other hand, the power spectrum receives 
contributions not only from the variance, but also from 
the covariance Ap(fe) Ap{— fc— e) between modes a small 
wavevector e apart. This covariance vanishes in the mean, 




Figure 3. Four-point configuration of wavevectors for the 
trispectrum in equation 1361 . which describes the covariance 
of power spectra of weighted densities. The short leg e, 
equation 1871 . produces a beat-coupling to large scales. 

but it couples to large-scale modes Ap(e) through quadratic 
nonlinearities. That is, the correlation between the product 
Ap{k)Ap{— fc— e) and the large-scale mode Ap(e) is the 
bispectrum 

(Ap(fc) Ap(-fc-e) Ap(e)) = B(k, -k-e, e) . (85) 

The bispectrum is zero for Gaussian fluctuations, but is 
driven away from zero by nonlinear gravitational growth. 

4.2 Tetrahedron 

The place where, prior to this section, we inadvertently 
discarded the large-scale beat-coupling, is equation 13711 . 
where we made the seemingly innocent approximation that 
the trispectrum T(fci, fc 2 , fc3, fc4) is a slowly varying function 
of what appears to be its arguments, fci to k±. This 
assumption is false, as we now show. 

For a statistically isotropic field (as considered in this 
paper), the trispectrum depends on six scalar arguments. 
This follows from the fact that a spatial configuration of 
four points is determined by the six lengths of the sides 
of the tetrahedron whose vertices are the four points. In 
Fourier space, the configuration is an object four of whose 
sides are equal to the wavevectors fci to k^. The object forms 
a closed tetrahedron (because ^ fci = 0), whose shape is 
determined by the six lengths of the sides of the tetrahedron. 

Figure |3] illustrates the configuration of interest in the 
present paper, that for the trispectrum in equation 13611 . 
Rewritten as a function of six scalar arguments, the 
trispectrum of equation 1361 is 

T(fci-fci,-fci-fci',fc 2 -fc 2 ,-fc 2 -fc 2 ) = (86) 
r(|fei— k[\, |fci+fc"|, |fc2— k' 2 \, I fc 2 -r- 2 1 , jfci— k 2 — fci— fc 2 |, e) 
where the wavevector e is defined by 
e = -(k' 1 +k'{)=k' 2 +k'i (87) 

which is small but not necessarily zero. The invalid 
approximation I37H is equivalent to approximating 

T(Jfci— fell , |fei+fe"|, \ k 2 — k' 2 \, | fc 2 —I— fc 2' 1 1 1 fci— fe2— fci— fei'li e ) 
« T(fci,fci,fc 2 ,fc 2 , |fci-fc 2 |,0) . (88) 
The problem with this approximation is apparent. Although 
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primed wavenumbers are small compared to unprimed ones, 
so that the approximation in the first five arguments 
is reasonable, in the last argument it is not valid to 
approximate a finite wavenumber e, however small, by zero. 
A valid approximation is, rather, 

T(\ki— fej |, |fei+fc"|, \k 2 — fc^l) |fca+fca l> \ki— k 2 — k'x— k 2 \, sj 
« T(k 1 ,k 1 ,k 2 ,k 2 ,\k 1 -k 2 \,e) . (89) 

As an example of the large-scale beat-coupling contribu- 
tions to the trispectrum that arise from the beat wavevector 
e, consider perturbation theory. 



4.3 Perturbation theory 

In perturbation theory (PT), the trispectrum can be split 
into snake and star contribut ions iScoccimarro et al]|l999t 
ISefnsatti fc RcoccimarrolEnnl 

T(fei,fe2,fe 3 ,fe 4 ) = 

4P(fci)P(fe)P(*i3)F a (fci, -*u)F a (fca, feu) 
+ cyclic (12 snake terms) 

+ P{k 1 )P{k 2 )P{k s ) [F 3 (fci, fc 2 , fcs) + perm. (6 terms)] 
+ cyclic (4 star terms) (90) 

where kij = ki + kj , and the second-order PT kernel F 2 is 
given by 



where v'i(k) is defined by equations 1491 and 1521 or 15311 . 

and the constant R a is 



„ ,. , . 5 x ( k\ k 2 \ 2 
F 2{kl ,k 2 ) = - + -(- + -) +- 



(91) 



with x = ki ■ k 2 . 

In the case of interest, where the trispectrum is that of 
equation 18611 . 4 of the 12 snake terms produce a coupling 
to large scales, those where the beat wavenumber fci3 in 
equation 1901 is small. In the (valid) approximation 1891 . 
the pertinent PT trispectrum is 

T[\ki— fer|, \ki+k" I, \k 2 —k' 2 \, \k 2 +k 2 \, \ki— k 2 — k[— k 2 |, e) 
w T(fci,fci,fc 2 ,fc2,|fe 1 -fe 2 |,0) 

+ 16 P(k 1 )P{k 2 )P(e)F 2 (k 1 ,-e)F 2 (k 2 ,e) (92) 

in which the term on the last line represents the large- 
scale beat-coupling contribution incorrectly ignored by the 
approximation 1881 . In equation 1361 for the covariance 
of shell-averaged power, this trispectrum, equation 1921 . is 
angle- averaged over the directions of fci and fc 2 . The angle- 
averaged second-order PT kernel is 



, , do k 17 
F2(fe ' £) 4^ = 21 



(93) 



and it follows that the last line of equation 1921 , when angle- 
averaged, is 16(17/21) 2 P(fei)P(fc 2 )P(e). 

Following the same arguments that led from equa- 
tion 1361 to equation 1401 . and then to equation 1541 . but 
with the beat-coupling term now correctly retained in the 
trispectrum, one finds that equation 1541 for the expected 
covariance of shell-averaged power spectra of weighted 
densities is modified to 



(A#(*i)A#(/Ba)) « fa {Ap(kx)Ap(k 2 )) 

+ 4 P a P(fc!)P(fc 2 ) «i( fe H (~ fe )F(fc) 



2.62 



(95) 



The reason for writing equation l|940 in this form, with the 
constant R a separated out, is that, as will be seen in H4.4I 
the same expression remains valid in the hierarchical model, 
but with R a the 4-point hierarchical snake amplitude. 

Figure [5] includes lines showing the predicted PT 
result for the variance of shell- averaged power of weighted 
density, equation 1361 . both with (solid lines) and without 
(dashed lines) beat-coupling. The PT variance with beat- 
coupling was obtained by numerically integrating the PT 
expression 1901 for the trispectrum 1861 in equation 1361 
(that is, without making the approximations 1891 or 1941 '). 
with the minimum variance weightings 1711 . and then 
multiplying by the factor 1//' = 2, equation 1741 . From 
this the PT variance without beat-coupling was obtained by 
setting P(e) = 0. The variance without beat-coupling agreed 
well with a direct PT evaluation of equation 1311 . 

Figure[5]shows that the beat-coupling contribution pre- 
dicted by perturbation theory seems to account reasonably 
well for the extra variance that appears at nonlinear scales 
in the weightings versus the ensemble method. 

We will return to equation 1941 in H4.5I below, but 
first consider the hierarchical model as a prototype of the 
trispectrum beyond perturbation theory. 

4.4 Hierarchical model 

Perturbation theory is valid only in the translinear regime. 
The behaviour of the trispectrum in the fully nonlinear 
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ICroton et a l. 2004) is consistent with a hierarchical model of 
highe r order correlations. In the hierarchical model jPeeblesI 
1980), the trispectrum is a sum of snake and star terms 

T(fe 1 ,fe 2 ,fe 3 ,fe 4 ) = 

Ra[P(k 1 )P(k 2 )P{ki 3 ) + cyclic (12 snake terms)] 

+ R b [P(ki)P(k 2 )P(k 3 ) + cyclic (4 star terms)] . (96) 

The PT trispectrum, equation 1901 . shows a hierarchical 
structure with hierarchical amplitudes R a and R D that 
are not constant, but rather depend on the shape of 
the trispectrum tetrahedron. At highly nonlinear scales, 
Scoccimarro & Frieman (1999) suggested an ansatz, dubbed 
hyperextended perturbation theory (HEPT), that the 
hierarchical amplitudes go over to the values predicted by 
perturbation theory for configurations collinear in Fourier 
space. For power law power spectra P(k) oc k n , HEPT 
predicts 4-point amplitudes 

54 - 27 2 n + 2 3" + 6 n 



Ra — Rb — 



(97) 



(94) 



2 (1 + 6 2" + 3 3 n + 6 6") 

As point e d ou t by IScoccimarro fc Friemanl II1999T) and 
iHamiltonl (I200CT) . HEPT is not entirely consistent because it 
predicts a covariance of power (Ap(fci) Ap(k 2 )) that violates 
the Schwarz inequality when > fe- 

In the hierarchical model with constant hierarchical 
amplitudes, 4 of the 12 snake terms produce a coupling to 
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large scales in the trispectrum of interest, equation 1861 . In 
the (valid) approximation 1891 , the hierarchical trispectrum 
is 

r(|fc 1 -fc' 1 |, |fei+fe'/|, |fc 2 -fc 2 |, |fca+fcal) Iki—ks-k'i-k'^e) 

» r(fci,fci,fc 2 ,/c2,|fei-fe 2 |,o) 

+ 4R a P(ki)P(k2)P(e) (98) 

in which the term on the last line represents the large-scale 
beat-coupling contribution. 

The hierarchical trispectrum 1981 looks similar to 
(slightly simpler than) the PT trispectrum 1921 . Following 
the same arguments as before, one recovers the same 
expression 1941 for the expected covariance of shell- averaged 
power spectra of weighted densities. 

4.5 Covariance of shell-averaged power spectra 
including large-scale coupling 

Suppose that either perturbation theory, H4.3I or the 
hierarchical model, H4.4I offers a reliable guide to the 
coupling of the nonlinear trispectrum to large scales, so 
that equation 1941 is a good approximation to the expected 
covariance of shell-averaged power spectra of weighted 
densities. 

Make the further assumption that the power spectrum 
is approximately constant over the large-scale wavevectors 
represented in v[ (fc) 

P(k) « P(2k b ) = constant for v-(k) ^ (99) 

where kb is the wavenumber at the box scale. The factor 
2 in 2kb in equation 11991 appears as a reminder that 
the wavevectors fc in i4(fc) are, equations 1491 and 1521 
or 1531 . sums of pairs of wavenumbers fc' represented in 
the weighting Wi(k'). For example, if the weightings are 
taken to be the minimum variance weightings given by 
equation 1711 . then kb = ki where ki is the wavenumber 
of the weighting. Approximation 1991 is in the same spirit 
as, but distinct from, the earlier approximation 1171 that 
the power spectrum is a slowly varying function. Note 
that equation 1991 does not require that P(0) ~ P(2kb) 
(which would certainly not be correct, because P(0) = 0), 
because ?4(0) is zero, which is true a priori in strategy one, 
equation 1521 , and ends up being true a posteriori in strategy 
two, equation 1531 . by the argument in H3.4I 

In the approximation 1991 . the summed expression on 
the right hand side of equation 1941 is 

^i4(fcX(-fc)p(fc) « f tJ P{2k b ) (ioo) 

k 

and equation 1941 reduces to 

(A#(fci)A#(fe,)) » fa [(Ap(fc!)Ap(fc 2 )) 

+ 4R a P(k 1 )P(k 2 )P(2k b )] (101) 

with the term on the last line being the large-scale beat- 
coupling contribution. 

Equation 11011 provides the fundamental justfication 
for the weightings method when beat-coupling is taken 
into account. It states that the covariance of shell- averaged 
power spectra of weighted densities is proportional to the 
sum of the true covariance (Ap(fci)Ap(fc2)) of shell- averaged 



power, and a beat-coupling term 4 R a P(ki)P(k2)P(2kb) 
proportional to power at (twice) the box wavenumber kb- 
The crucial feature of equation 11011 is that the constant 
of proportionality fa, equation 1551 . depends only on the 
weightings Wi(k), and is independent either of the power 
spectrum P(k) or of the wavenumbers ki and k%. 

In the limit of infinite box size, the beat-coupling 
contribution to the covariance of power spectra of weighted 
densities in equation 11011 goes to zero, P(2kb) — > P(0) = 
as kb — > 0, and the covariance becomes proportional to 
the true covariance (Ap(fci)A^(fe)) of power. However, in 
cosmologically realistic simulations, such as illustrated in 
Figure [5] and discussed further in |SJ the beat-coupling 
contribution, far from being small, is liable to dominate at 
nonlinear scales. 

Beyond perturbation theory or the hierarchical model, 
the weightings method remains applicable just so long as the 
hierarchical amplitude R a in equation 11011 is independent 
of the weightings ij. In general, R a could be any arbitrary 
function of ki, ki, and the box wavenumber kb- 

4.6 Not quite minimum variance weightings 

Section |3 derived sets of minimum variance weightings valid 
when the covariance, and the covariance of covariance, of 
power spectra of weighted densities took the separable forms 
given by equations 1541 and 1591 . When beat-scale coupling 
is included, the covariance of power, equation 11011 . still 
takes the desired separable form (as long as the hierarchical 
amplitude R a is independent of the weightings ij), but the 
covariance of covariance of power (eq. 1A31 of Appendix [XJ 
does not. 

In Appendix ^ we discuss what happens to the 
minimum variance derivation of J|] when beat-coupling is 
included. We argue that the minimum variance weightings of 
>!3.3l are no longer exactly minimum variance, but probably 
remain near minimum variance, and therefore fine to use in 
practice. 

The factor 2 on the right hand side of equation 1771 is 
no longer correct when beat-coupling is included, but may 
remain a reasonable approximation. 



5 DISCUSSION 

As shown in [Jlj the covariance of nonlinear power receives 
beat-coupling contributions from large scales whenever 
power is measured from Fourier modes p(k) that have 
a finite spread in wavevector fc, as opposed to being 
delta-functions at single discrete wavevectors. Physically, 
the large-scale beat-coupling arises because a product 
Ap(k) Ap(— fc— e) of Fourier amplitudes of closely spaced 
wavevectors couples by nonlinear gravitational growth to the 
beat mode Ap(e) between them. 

The beat-coupling contribution does not appear when 
covariance of power is measured from ensembles of periodic 
box simulations, because in that case power is measured 
from products of Fourier amplitudes Ap(fc) Ap(— fc) at single 
discrete wavevectors. Here the "beat" mode is the mean 
mode, fc— fc = 0, whose fluctuation is by definition always 
zero, Ap(0) = 0. 

There is on the other hand a beat-coupling contribution 
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when covariance of power is measured by the weightings 
method, because the Fourier modes of weighted density are 
spread over more than one wavevector. 

For weightings constructed from combinations of fun- 
damental modes, as recommended in JJ] the covariance of 
power spectra of weighted densities receives beat-coupling 
contributions from power near the box fundamental kb- 
The beat-coupling and normal contributions to the variance 
i(fc) 2 ) of nonlinear power are in roughly the ratio 
P(2kb)/P(k) of power at the box scale to power at the 
nonlinear scale, according to equation 11011 . 

In cosmologically realistic simulations, box sizes are 
typically around the range 10 2 -10 3 ft _1 Mpc. This is just the 
scale at which the power spectrum goes through a broad 
maximum. For example, in observationally concordant 
ACDM models, power goes thro ugh a broad maximum 
at fc pcak w 0.016/iMpc" 1 (e.g. iTeemark et al.1 l2004al) . 
corresponding to a box size 47r/fc poa k ~ 800 h~ Mpc. Power 
at the maximum is about 25 times greater than power at 
the onset fctrans ~ 0.3 h Mpc -1 of the translinear regime, 
-P(fc P eak)/^ > (fctr a n S ) « 25, and the ratio P(fc poak ) / P{k) of 
power increases at more nonlinear wavenumbers k. 

It follows that in cosmologically realistic simulations 
the beat-coupling contribution to the covariance of power 
is liable to dominate the normal contribution. This is con- 
sistent with the numerical results illustrated in Figure Hand 
discussed by iRimes fc Hamilton! <|2006T) . which show that 
the variance of power measured by the weightings method 
(which includes beat-coupling contributions) substantially 
exceeds, at nonlinear scales, the variance of power measured 
by the ensemble method (which does not include beat- 
coupling contributions). 

5.1 Relevance to real galaxy surveys 

In real galaxy surveys, measured Fourier modes inevitably 
have finite width [Afe| ~ 1/R, where R is a characteristic 
linear size of the survey. The characteristic size R varies from 
100 ft Mpc to a few 1000 /i _1 Mpc (an upper limit is set by 
the comoving horizon distance, which is about 10 4 /i _1 Mpc 
in the concordant ACDM model). 

It follows that the covariance of nonlinear power 
measured in real galaxy surveys is liable to be dominated not 
by the "true" covariance of power (the covariance of power 
in a perfect, infinite survey), but rather by the contribution 
from beat-coupling to power at the scale of the survey. 

This means that one must take great care in using nu- 
merical simulations to estimate or to predict the covariance 
of nonlinear power expected in a galaxy survey. The scatter 
in power over an ensemble of periodic box simulations 
will certainly underestimate the covariance of power by a 
substantial factor at nonlinear scales, because of the neglect 
of beat-coupling contributions. 

A common and in principle reliable procedure is to 
estimate the covariance of power of a galaxy survey from 
mock surveys "observed" with the same selection rules as 
the real survey from numerical simulati ons large enough 
to encompass the en ti re survey (e.g . |C 1 oiletalJ_|20 i 0jt 
Padilla fc Baughl 120031: Teemark et*aH l2004al: IVan et all 
20041: iBlaizot et alJl20Q.4 iFrith et alJl2005l : ICole et al.ll200,4 
Eisenstein et alJl200.4 IPark et al.ll2005h . 

It is important that numerical simulations be genuinely 



large enough to contain a mock survey. One should be 
wary about estimating covariance of power from mock 
surveys ext racted from small p eriodic boxes replicated many 
times (e.g. lYang et al.ll2004Ti . since such boxes are liable 
to be missing power at precisely those wavenumbers, the 
inverse scale size of the mock survey, where beat-coupling 
should in reality be strongest. Beat-coupling arises from a 
real gravitational coupling to large scale modes, and the 
simulation from which a mock survey is extracted must be 
large enough to contain such modes. 

Further, it would be wrong to take, say, a volume- 
limited subsample of a galaxy survey, and then to estimate 
the covariance of power from an ensemble of periodic nu- 
merical simulations whose size is that of the volume-limited 
subsample. A volume- limited subsample of observational 
data retains beat-coupling contributions to the covariance 
of power, whereas periodic box simulations do not. 



6 SUMMARY 

This paper falls into two parts. In the first part, JUand J|] 
we proposed a new method, the weightings method, that 
yields an estimate of the covariance of the power spectrum 
of a statistically homogeneous and isotropic density field 
from a single periodic box simulation. The procedure is 
to apply a set of weightings to the density field, and to 
measure the covariance of power from the scatter in power 
over the ensemble of weightings. In J2| we developed the 
formal mathematical apparatus that justifies the weightings 
method, and in we derived sets of weightings that 
achieve minimum variances estimates of covariance of power. 
Section 13.61 gives a step-by-step recipe for applying the 
weightings method. We recommend a specific set of 52 
minimum variance weightings containing only combinations 
of fundamental modes. 

In the second part of this paper, 21 and |5] we discuss 
an unexpected glitch in the procedure, that emerged from 
the periodic box numerical simulations d escribed in the 
companion paper (IRimes fc Hamiltonll2006l) . The numerical 
simulations showed that, at nonlinear scales, the covariance 
of power measured by the weightings method substantially 
exceeded that measured over an ensemble of independent 
simulations. 

In [21 we argue from perturbation theory that the 
discrepancy between the weightings and ensemble methods 
arises from "beat-coupling", in which products of closely 
spaced Fourier modes couple by nonlinear gravitational 
growth to the large-scale beat mode between them. Beat- 
coupling is present whenever nonlinear power is measured 
from Fourier modes that have a finite spread in wavevector, 
as opposed to being delta-functions at single discrete 
wavevectors. Beat-coupling affects the weightings method, 
because Fourier modes of weighted densities have a finite 
width, but not the ensemble method, because the Fourier 
modes of a periodic box are delta-functions of wavevector. 

As discussed in beat-coupling inevitably affects real 
galaxy surveys, whose Fourier modes necessarily have a 
finite width of the order of the inverse scale size of the survey. 
Surprisingly, at nonlinear scales, beat-coupling is liable to 
dominate the covariance of power of a real survey. One would 
have thought that the covariance of power at nonlinear scales 
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would be dominated by structure at small scales, but this is 
not true. Rather, the covariance of nonlinear power is liable 
to be dominated by beat-coupling to power at the largest 
scales of the survey. 

A common and valid procedure for estimating the 
covariance of power from a real survey is the mock survey 
method, in which artificial surveys are "observed" from 
large numerical simulations, with the same selection rules 
as the real survey. It is important that mock surveys be 
extracted from genuinely large simulations, not from many 
small periodic simulations stacked together, since stacked 
simulations miss the large-scale power essential to beat- 
coupling. 

Finally, it should be remarked that, although this paper 
has considered only the covariance of the power spectrum, 
it is likely that, in real galaxy surveys and cosmologically 
realistic simulations, beat-coupling contributions dominate 
the nonlinear variance and covariance of most other 
statistical measures, including higher order n-point spectra 
such as the bispectrum and trispectrum, and n-point 
correlation functions in real space, including the 2-point 
correlation function. 
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Figure Al. Eight-point configuration of wavevectors contribut- 
ing to the covariance of covariance of power spectra of weighted 
densities, equation 1A11 . The central tetrahedron of short legs 
produces beat-couplings to large scales. 



APPENDIX A: MINIMUM VARIANCE 
WEIGHTINGS AND BEAT-COUPLING 

This Appendix describes how the beat-coupling contribu- 
tions to covariance of power discussed in &\ modify the 
minimum variance arguments in Sj3| The conclusion is that 
the minimum variance weightings given in HH.HI are no longer 
exact minimum variance, but are probably near minimum 
variance, and therefore fine to use in practice. 

With no approximations at all, the covariance of 
covariance of shell-averaged power spectra of weighted 
densities (with the deviations Ap^(fc) in power being taken 
relative to the measured rather than the expected mean 
power, eqs. 1471 or 1481 1. takes the generic form 

[A#(*i)A#(fca) - (A#( fcl )A#(fca))] 

x [A#(fcs)A#(*4) - (A#(fcs)A#(fc4)); 
= ^ Ui(ei)Wj(e2)«j-(e3)«K e 4) (Al) 

ei-t-e2+e3+£4— ° 

£/(fcl — fel , — fel — fel , k^. — ^4; — k^ — ^4 ) 



E 



where v'^k) is defined by equations l|49|l and (l"i2> or l)53[l . 

the wavevectors e„ are defined by 



= —k'—k'L for n = 1, 



(A2) 



and E is an 8-point object, a sum of products of n- 
point functions adding up to 8 points, as enumerated in 
equation (1431 1 . 

Figure IA1I illustrates the configuration of the 8-point 
function that contributes to the 8-point object E in 
equation HA It . The short legs ei, £4 of the configuration 
constitute a tetrahedron, whose 6 sides generate beat- 
couplings to large scale. None of the legs e„ is zero, because 
a zero leg would make zero contribution to equation 1AH . 



since t^(0) = 0, equation (I52H or, a posteriori, equation (15311 : 
but it is possible for the sum of a pair of the short legs to 
be zero. It is these configurations, where the sum of a pair 
of short legs is zero, that prevent equation 1A11 from being 
separated, as in equation 1591 . into a product of a factor 
g'ij that depends only on the weightings and a factor that 
is independent of the weightings. Note that it is only the 8- 
point function itself, not lower-order functions, that prevent 
separability: lower-order functions depend on at most three 
of the four e n , and a triangle with three non-zero sides has 
(of course) no zero sides. 

In either perturbation theory or the hierarchical model, 
and in the various valid approximations made in this paper 
(firstly, that n-point spectra are slowly varying functions 
of their arguments except that small arguments £ are not 
replaced by zero, and secondly, that shells are broad), the 
covariance of covariance of shell- averaged power spectra of 
weighted densities, equation HAIL reduces to 



[A^(fcx)A^(*a) - (A^(*i)A^(*a)>] 

x [A^(fc 3 )A^(fc 4 ) - (/\p' ] {k 3 )Ap' J {k A ))\ 



^9ij ~ f'fiifjj — v fi; 



(A3) 



where fL and g'^ are given by equations 155H and (16011 . The 
quantities A, fi, v in equation <A3t are each functions of 
fcij &2, &3, &4, and the box wavenumber kb, but, importantly, 
are independent of the weightings ij. The f'ufjj (respectively 
f[f) term in equation <A3t arises from terms where 
£i+£2 = (respectively £i+£3 = or £i+£4 = 0) 
in equation HAIL All three terms on the right hand 
side of equation <A3t contain large-scale beat-coupling 
contributions, proportional to one, two, or three factors of 
large-scale power. The second (ji) and third (y) terms in 
equation JAlfc are written with negative signs because their 
effect is such as to cancel some of the beat-coupling terms 
appearing in the first (A) term (that is, some of the beat- 
coupling terms proportional to P(2kb) in A should really be 
proportional to P(0) = 0; the p and v terms remove these 
terms). It is to be expected that A, A — fi, and A — u are all 
positive. At linear scales, where fluctuations are Gaussian, 
the beat-couplings generated by nonlinear evolution are 
small, so that A ^> fi, v. At nonlinear scales, however, [i 
and v could be an appreciable fraction of A. 

The derivation of minimum variance weightings in Sj^l 
involved summing over weighings ij, equation l|68^ . Consider 
the corresponding double sum over ij of equation (|A3(l . The 
sum over yields the same result as before, equation 169H . 
Adjoining the sum over f'ufjj from equation l|A3fl modifies 
equation I69t to 



TaK2 E ~ (V/X)flifjj] 



(f'N) 



1-3 
1 



n'(0) 2 



(i-£)t/(o) 2 + 5>'(fc)| 



fe^O 



(A4) 



because f' u = /' = u'(0) for all i. The minimum variance 
weightings given in H3.3l were absolute minimum variance in 
the sense that u'(fe) = for k 7^ 0. The same minimum 
variance weightings continue to achieve absolute minimum 



variance for equation HA4I . reducing its right hand side to 
the irreducible minimum 1 — fi/X. 

Thus the minimum variance weightings of H3. 31 remain 
minimum variance as long as only the first two terms (A 
and fi) of equation 1A3II are considered. The third (v) term 
breaks the minimum variance derivation. However, this third 
term is likely to be subdominant compared to the first 
two. The quantity fL in the third term of equation 1A3I is 
proportional to the covariance of power between weightings i 
and j, equation I10H , and the Schwarz inequality guarantees 
that 

fij ^ fafjj (A5) 

so that there is a natural tendency for the third term of 
equation IA3t to be dominated by the second. The only 
way for the third term to be large is for the power spectra 
from different weightings ij to be highly correlated with each 
other. Physically, however, the most accurate estimate of 
covariance of power should come from averaging over many 
uncorrelated weightings, in which case fl? <C fLf'jj f° r most 
weightings i 7^ j. Thus, as just stated, it is to be expected 
that the third term should be subdominant compared to the 
first two. 

In summary, to the extent that either perturbation 
theory or the hierarchical model provide a reliable guide to 
the behaviour of high-order correlations, and to the extent 
that the third term of equation IA3t is subdominant, as it 
should be, the minimum variance weightings of >!3.3l should 
remain near minimum variance, good enough for practical 
application. 
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